%% plotting the figures for the PP paper

%loading the spectra data


load /nfs/satmag_work/mnair/projects/longp/alldays JULI_SEG ACE_SEG;
JULI_SEG = JULI_SEG.*24.366*1e-3; %mV/m
[Txy_short,F_short] = tfestimate(ACE_SEG,JULI_SEG,hanning(72),0,72,1/(5*60)); %Txy is agains calculated just to get
[Cxy_short,F_short] = mscohere(ACE_SEG,JULI_SEG,hanning(72),0,72,1/(5*60));
N_data = length(ACE_SEG) / 72;
Err_short = sqrt( 1/(2*(N_data-1)) .* ( (1-Cxy_short)./Cxy_short ) ) .* abs(Txy_short);
Err_short_lg = 0.434 * Err_short ./ abs (Txy_short);

%% Long period data

load /nfs/satmag_work/mnair/projects/longp/champ_eef_TF_data options ACE_SEG ACE_SEG_EZ TIME_SEG CHAMP_SEG N_data N_seg F_long Txy_long Cxy_long Err_long;
[Cxy_long,F_long] = mscohere(ACE_SEG, CHAMP_SEG,hanning(options.len),0,options.len,1/(options.des_int*3600)); %1/(5*60) = sampling frequency in Hz )
[Pxy_long,F] = cpsd(CHAMP_SEG,ACE_SEG,hanning(options.len),0,options.len,1/(options.des_int*3600)); %
[Pxx_long,F] = pwelch(CHAMP_SEG,hanning(options.len),0,options.len,1/(options.des_int*3600));
[Txy_long,F_long] = tfestimate(ACE_SEG,CHAMP_SEG, hanning(options.len),0,options.len,1/(options.des_int*3600));

%Compute the TF magnitude error
Err_long = sqrt( 1/(2*(N_data-1)) .* ( (1-Cxy_long)./Cxy_long ) ) .* abs(Txy_long);
Err_long_lg = 0.434 * Err_long ./ abs (Txy_long);

% TF magnetitude EEF and IEF Ez
[Tez_long,F_long] = tfestimate(ACE_SEG_EZ,CHAMP_SEG, hanning(options.len),0,options.len,1/(options.des_int*3600));
% Coherence between EEF and IEF Ez
[Cez_long,F_long] = mscohere(ACE_SEG_EZ,CHAMP_SEG, hanning(options.len),0,options.len,1/(options.des_int*3600));
%Error
Err_Ez_long = sqrt( 1/(2*(N_data-1)) .* ( (1-Cez_long)./Cez_long ) ) .* abs(Tez_long);
Err_Ez_long_lg = 0.434 * Err_Ez_long ./ abs (Err_Ez_long);


%errorbar(log10((1./F_long(2:end))./(3600)), log10(abs(Txy_long(2:end))),Err_long_lg(2:end)/2,Err_long_lg(2:end)/2,'c');

%% 
%%Coherence

load  /nfs/satmag_work/mnair/projects/longp/JULIA_GT20_FINAL;


figure1 = figure;
axes('Parent',figure1,'XTick',[0.1 .7 10 40],...
    'XScale','log',...
    'XMinorTick','on');
set(gca,'FontSize',16);
box('on');
hold('all');
xlabel('Period');
ylabel('coherence');
hold on;
%confidence level of coherence
%ref Thompson, R.O., 1979: Coherence Significance Levels. J. Atmos. Sci.,
%*36*, 2020�2021.
alpha = 0.95;
ci = 1- (1-alpha)^(1/((N_data)-1));
aa=axis;
h=line([aa(1) aa(2)],[ci,ci],'LineStyle','-.','color','b');
signif_coh = find(Cxy_short(2:end) > ci);
semilogx((1./(3600*F_short(signif_coh))),Cxy_short(signif_coh),'b','LineWidth',2);
axis([0.1,40,0,0.8]);
hold on;

load /nfs/satmag_work/mnair/projects/longp/CHAMP_LONG_PERIOD_GT20_FINAL

semilogx((1./(3600*F_long(2:end))),Cxy_long(2:end),'r','LineWidth',2);

%set(gcf,'position',[1002 484 784 583]);
set(gca,'XTick',[6/60,10/60,20/60,30/60,1,2,4,6,10, 30]);
set(gca,'XTickLabels',[' 6';'10';'20';'30';' 1';' 2';' 4';' 6';'10';'30']);
%text(6/60,-0.08,'|<-','FontSize',16);
%text(8,-0.08,'->|','FontSize',16);
%text(40/60,-0.08,'-> | <-','FontSize',16);
text(10/60,-0.08,'minutes','FontSize',16);
text(10,-0.08,'hours','FontSize',16);

alpha = 0.95;
ci = 1- (1-alpha)^(1/((N_data)-1));
aa=axis;
h=line([6,40],[ci,ci],'LineStyle','-.','color','r');



%% tf
figure1=figure;
axes('Parent',figure1,'XTick',[0.1 1 40],...
    'XScale','log',...
    'XMinorTick','on');
set(gca,'FontSize',16);
box('on');
hold('all');
xlabel('Period');
ylabel('Magnitude (dB)');
hold on;

Xmag_short = abs(Txy_short);          % Spectral magnitude
Xdb_short = 20*log10(Xmag_short);   % Spectral magnitude in dB

Xmag_long = abs(Txy_long);          % Spectral magnitude
Xdb_long = 20*log10(Xmag_long);   % Spectral magnitude in dB

Xmag_long_ez = abs(Tez_long);          % Spectral magnitude
Xdb_long_ez = 20*log10(Xmag_long_ez);   % Spectral magnitude in dB

% XdbMax = max(Xdb);      % Peak dB magnitude
% Xdbn = Xdb - XdbMax;    % Normalize to 0dB peak
%
% dBmin = -100;           % Don't show anything lower than this
% Xdbp = max(Xdbn,dBmin); % Normalized, clipped, dB mag spec

L = Cxy_long > options.coh_plot_lim & (1./F_long)./(3600) >= options.minimum_period_plot ;

semilogx((1./(3600*F_short(2:end))),Xdb_short(2:end),'bs-','LineWidth',2);
hold on;
semilogx((1./(3600*F_long(L))),Xdb_long(L),'rs-','LineWidth',2);

axis([0.1 40 -55 -20 ]);
%set(gcf,'position',[1002 484 784 583]);
set(gca,'XTick',[6/60,10/60,20/60,30/60,1,2,4,6,10, 30]);
set(gca,'XTickLabels',[' 6';'10';'20';'30';' 1';' 2';' 4';' 6';'10';'30']);
%text(6/60,-0.08,'|<-','FontSize',16);
%text(8,-0.08,'->|','FontSize',16);
%text(40/60,-0.08,'-> | <-','FontSize',16);
aa=axis;
text(aa(1),aa(3),'minutes','FontSize',16);
text(aa(2),aa(3),'hours','FontSize',16);


h=line([(1./(3600*F_short(2:end))) (1./(3600*F_short(2:end)))]',[Xdb_short(2:end) ...
    - 20*Err_short_lg(2:end)  Xdb_short(2:end) + 20*Err_short_lg(2:end) ]','color','b','LineWidth',2)
h=line([(1./(3600*F_long(L))) (1./(3600*F_long(L)))]',[Xdb_long(L) ... 
    - 20*Err_long_lg(L)  Xdb_long(L) + 20*Err_long_lg(L) ]','color','r','LineWidth',2);


%% Create final filter coefficients from JULIA and CHAMP data
% load /nfs/satmag_work/mnair/projects/longp/julia_champ_ace_tf_data_final ... 
%     Txy_short Cxy_short Txy_long Cxy_long F_short F_long Err_long_lg Err_short_lg Err_short Err_long;

load  /nfs/satmag_work/mnair/projects/longp/JULIA_GT20_FINAL ; % Suppoesed to be final :-)
load /nfs/satmag_work/mnair/projects/longp/CHAMP_LONG_PERIOD_GT20_FINAL


L = 1./(F_short*3600) >= 5.9999  | F_short == 0 | 1./(F_short*3600) <=    (10/60) ; % eleminate unwanted frequencies

% L = ones(size(F_short)); % select only a few frequencies
% L([3,5,8,13,18,23,28,40]) = 0;
% L = logical(L);

Txy_short(L) = []; % remove the unselected frequencies
Cxy_short(L) = [];
F_short(L) = [];
Err_short_lg(L) = [];
Err_short(L) = [];

% remove periods less than 6 hours and DC component from CHAMP data

L = 1./(F_long*3600) < 6  | 1./(F_long*3600) > 40 ;
%L = 1./(F_long*3600) < 6 ;

Txy_long(L)  = [];% remove the unselected frequencies
Cxy_long(L) = [];
F_long(L) = [];
Err_long_lg(L) = [];
Err_long(L) = [];

[F_all,IA,IB] = union(F_long,F_short); % merge JULIA and CHAMP magnitudes
[A,B] = sort([F_long(IA);F_short(IB)]);

%A and F_all should be the same;

T_all = [Txy_long(IA);Txy_short(IB)]; % Final TF
C_all = [Cxy_long(IA);Cxy_short(IB)]; %Final coherence
Error_sn = [Err_long_lg(IA);Err_short_lg(IB)]; %S/N
Error_mag = [Err_long(IA);Err_short(IB)]; %TF Error magnitude
T_all = T_all(B); %sort the Txy according the the frequency mix;
C_all = C_all(B); 
Error_sn = Error_sn(B);
Error_mag = Error_mag(B);
F_scaled = pi*(F_all./max(F_all)); %scale the frequency between o and pi
%T_all(1) = 0;
%  Error_sn(2:3) = Error_sn(2:3)/100;
%  Error_sn(6:10) = Error_sn(8:12)/100;
%  %Error_sn(20:25) = Error_sn(20:25)/50;
%  Error_sn(35:38) = Error_sn(35:38)/100;
C_all(2:3) = 5.0 ;
C_all(6:7) = 5.0;
C_all(30:end) = 2;



%%
close all;
 open '/nfs/satmag_work/mnair/projects/longp/Figures/transfer_function_new.fig';
 h211 = gca;
 hold on;
warning off;

sa = 3;

for jj = 1:sa     % maximum number of poles
    for ii =1:jj  % number of zeros is always less than number of poles
               
       
%[b_n,a_n]=invfreqz(T_all,F_scaled,ii,jj,[],300,0.01); %No weight
[b_n,a_n]=invfreqz(T_all,F_scaled,ii,jj,C_all,300,0.001);
%[b_n,a_n]=invfreqz(T_all,F_scaled,ii,jj,1./Error_sn, 300, 0.001); 

% Filter fit
[h_n] = freqz(b_n,a_n,F_all,1/(5*60));

%find the maximum difference between observed tf and filter ts
maxi=max(abs(db(abs(h_n(2:end)))-db(abs(T_all(2:end)))));

% The floowing 2 lines uses the analog filter fit, which didn't work well
% [b_n,a_n]=invfreqs(T_all,F_all,4,2,Error_sn,3000,0.0001);
% [h_n] = freqs(b_n,a_n,F_all);
% subplot(211);
% loglog((1./(3600*F_all)),abs(T_all),'rs','MarkerSize',5);
% hold on;
axes(h211);
h = loglog((1./(3600*F_all)),abs(h_n),'b','LineWidth',2);
% set(gca,'FontSize',16);

% for i = 1: length(T_all);
%     plot([1./(3600*F_all(i)) 1./(3600*F_all(i)) ], [abs(T_all(i)) - ... 
%         Error_mag(i) abs(T_all(i)) + Error_mag(i)],'r','LineWidth',1);
%     hold on;
% end;
%     
% axis([0.1 40 1e-3 1e-1 ]);

set(gca,'XTick',[6/60,10/60,20/60,30/60,1,2,4,6,10, 30]);
set(gca,'XTickLabels',[' 6';'10';'20';'30';' 1';' 2';' 4';' 6';'10';'30']);

% aa=axis;
% text(aa(1),aa(3),'minutes','FontSize',16);
% text(aa(2),aa(3),'hours','FontSize',16);
% title ( sprintf('%d-%d',ii,jj));
% saveas(gcf,['/data/backup/mnair/temp/fliterresp' sprintf('%d-%d',ii,jj)],'png');
rms1 = sqrt(mean((abs(T_all)-abs(h_n)).^2));
title ( sprintf('%d-%d error(db) = %6.4f',ii,jj, maxi));

len = 120;
h212 = subplot(212);
step = zeros([1,len]);
step(13:end) = 1; % positive step
hh=plot((1:5:len*5)./60,filter(b_n,a_n,step).*1.0,'r-','LineWidth',1);
hold on;
plot((1:5:len*5)./60,ppefm(step,'Ey').*1.0,'b-','LineWidth',1);
rms2 = sqrt(mean(filter(b_n,a_n,step)-ppefm(step,'Ey')).^2);
title ( sprintf('rms = %6.4f',rms2));
grid on;

axis([0 10 0 0.07]);
fprintf('%d %d rms1 %f rms2 %f\n',ii,jj,rms1,rms2);
%if rms1 < 0.01 & rms2 < 0.001,
saveas(gcf,['/data/backup/mnair/temp/eefresp' sprintf('%d-%d',ii,jj)],'png');
%end;
%close all
delete(h);
delete(hh);
    end;
end;

%% Plot responses to synthetic IEF functions

%a=a_n;b=b_n;

         a = [ 1.0000000000 -0.5857685663 -0.9552092349 0.6300518820];
         b = [  0.0179652952 0.0092006195 -0.0177488188 -0.0089805341];

        ief_amp = 1; % IEF amplitude in mV/m
len = 120;
step = zeros([1,len]);
step(13:36) = ief_amp; % positive  box
figure(2);
plot((1:5:len*5)./60,step,'b-','LineWidth',4);
grid on;
hold on;
figure(3);
plot((1:5:len*5)./60,filter(b,a,step).*1.0,'b.-','LineWidth',4);%1.5 was for tf made only from Julia
hold on;
grid on;


step = zeros([1,len]);
step(13:end) = ief_amp; % positive step
figure(2);
plot((1:5:len*5)./60,step,'r-','LineWidth',2);
figure(3);
plot((1:5:len*5)./60,filter(b,a,step).*1.0,'r-','LineWidth',2)

step = zeros([1,len]);
step(25) = ief_amp; % dac delta
figure(2);
plot((1:5:len*5)./60,step,'k-','LineWidth',3);
hold on;
figure(3);
plot((1:5:len*5)./60,filter(b,a,step).*1.0,'k-','LineWidth',3)
grid on;

step = zeros([1,len]);%
step(30:40) = triang(11)'*ief_amp;
figure(2);
plot((1:5:len*5)./60,step,'c-','LineWidth',2);
hold on;
figure(3);
plot((1:5:len*5)./60,filter(b,a,step).*1.0,'c-','LineWidth',2)
grid on;


figure(2);
set(gca,'FontSize',16);
axis([0,10,-1,2]);
xlabel('Time (hours)')
ylabel('IEF Ey mV/m');


figure(3);
set(gca,'FontSize',16);
%axis([0,6,-0.1,0.1]);
axis([0,10,-0.1,0.1]);
%set(gca,'Ytick',[-0.05, -0.025,0,0.025,0.05]);
xlabel('Time (hours)');
ylabel('Equatorial Zonal EF mV/m');

%% Plot Filter-fit

 [h,w] = freqz(b,a,F_all,1/(5*60));
 [h_n] = freqz(b_n,a_n,F_all,1/(5*60));
        
        subplot(212);
%         semilogx(F_all,real(h),'b',F_all,imag(h),'b');hold on; 
%         semilogx(F_all,real(h_n),'k',F_all,imag(h_n),'k')    ;
%         semilogx(F_all,real(T_all),'r.',F_all,imag(T_all),'ro');
%         semilogx(F_all(L),real(T_all(L)),'r*',F_all(L),imag(T_all(L)),'r*');

        loglog(F_all,abs(h),'b');hold on; 
        loglog(F_all,abs(h_n),'k');   ;
        loglog(F_all,abs(T_all),'r.');
        loglog(F_all(L),abs(T_all(L)),'r*');

%% plot ISR-ACE correlation
load /nfs/satmag_work/mnair/projects/longp/data_for_isr_ace_correlation_plot xplot corplot dataplot corerr

errorbar(xplot, corplot./dataplot,corerr./dataplot, corerr./dataplot ,'ks');


xticks = 0:2:24 ;

set(gca,'XTick', xticks);


xticks_lt = xticks - 5;

xticks_lt (xticks_lt <= 0) = xticks_lt (xticks_lt <= 0) + 24;

set(gca,'XTickLabel',  reshape(sprintf('%2.0f',xticks_lt),[2,length(xticks_lt)])' );

set(gca,'FontSize',16);
xlabel('Local Time (Hours)');
ylabel('Correlation');
legend('Jicamarca ISR - ACE IEF Ey');
grid on;
hold on

% plot the polyfit

xplot_polyfit = [-fliplr(xplot)-0.5 xplot xplot+24];

%plot([-fliplr(xplot)-0.5 xplot xplot+24] ,[corplot corplot corplot ]./[dataplot dataplot dataplot ])

[p,s ] = polyfit(xplot_polyfit ,[corplot corplot corplot ] ./ [dataplot dataplot dataplot ], 21);


% plot(xplot_polyfit ,...
%     [corplot(end-2:1:end) corplot corplot(1:3) ]./ ...
%     [dataplot(end-2:1:end) dataplot dataplot(1:3) ]);


%
y = polyval(p,xplot);
plot(xplot, y,'r','LineWidth',2);
axis([-1,25,-.8,.8]);

%% plot LT scaling

load /nfs/satmag_work/mnair/projects/longp/data_for_isr_ace_correlation_plot xplot corplot dataplot corerr

factor = 1/max(y);
correlation_data = [corplot./dataplot corplot./dataplot corplot./dataplot] *  factor ;
[p,s ] = polyfit(local_t_xscale ,correlation_data, 21);
y = polyval(p,0:1:24);
axis([-1,25,-1,1]);

grid on
set(gca,'FontSize',16);
xlabel('Local Time (Hours)');
ylabel('LT factor')


%% plot IEF Ez transfer function magnitude

   
load '/data/backup/mnair/ace_tensor/Electric_Field_Matrix' seg_length EEF_SEG IEF_EY_SEG IEF_EZ_SEG TIME_SEG;
EEF_SEG = EEF_SEG.*24.366*1e-3; %mV/m
[Cxx_IEF_Ez,F] = mscohere(EEF_SEG,IEF_EZ_SEG,hanning(seg_length),0,seg_length,1/(5*60)); %1/(5*60) = sampling frequency in Hz (5*60 = 300 seconds)
[Ty,F_short] = tfestimate(IEF_EZ_SEG,EEF_SEG,hanning(seg_length),0,seg_length,1/(5*60));
N_data = length(IEF_EZ_SEG)/ seg_length;
Err_Ey_short = sqrt( 1/(2*(N_data-1)) .* ( (1-Cxx_IEF_Ez)./Cxx_IEF_Ez ) ) .* abs(Ty);
Err_Ey_short_lg = 0.434 * Err_Ey_short ./ abs (Ty);

%load  /nfs/satmag_work/mnair/projects/longp/champ_eef_tf_apgt20 CHAMP_SEG ACE_SEG N_data len des_int;
load /nfs/satmag_work/mnair/projects/longp/champ_eef_TF_data options ACE_SEG ACE_SEG_EZ TIME_SEG CHAMP_SEG N_data N_seg F_long Txy_long Cxy_long Err_long;
% TF magnetitude EEF and IEF Ez
[Tez_long,F_long] = tfestimate(ACE_SEG_EZ,CHAMP_SEG, hanning(options.len),0,options.len,1/(options.des_int*3600));
% Coherence between EEF and IEF Ez
[Cez_long,F_long] = mscohere(ACE_SEG_EZ,CHAMP_SEG, hanning(options.len),0,options.len,1/(options.des_int*3600));
%Error
Err_Ez_long = sqrt( 1/(2*(N_data-1)) .* ( (1-Cez_long)./Cez_long ) ) .* abs(Tez_long);
Err_Ez_long_lg = 0.434 * Err_Ez_long ./ abs (Err_Ez_long);


open '/nfs/satmag_work/mnair/projects/longp/Figures/transfer_function.fig';
% open the IEF Ey figures and plot the IEF Ez on to it. Do not save this
% file in the same name. 
hold on;
Xmag1 = abs(Tez_long);          % Spectral magnitude
Xmag2 = abs(Ty);
Xdb1 = 20*log10(Xmag1);   % Spectral magnitude in dB
Xdb2 = 20*log10(Xmag2);

semilogx((1./(3600*F_long(2:end))),Xdb1(2:end),'rs-','LineWidth',2);
semilogx((1./(3600*F(2:end))),Xdb2(2:end),'bs-','LineWidth',2);

h=line([(1./(3600*F_short(2:end))) (1./(3600*F_short(2:end)))]',[Xdb2(2:end) ...
    - 20* Err_Ey_short_lg(2:end)  Xdb2(2:end) + 20* Err_Ey_short_lg(2:end) ]','color','b','LineWidth',2);

h=line([(1./(3600*F_long(2:end))) (1./(3600*F_long(2:end)))]',[Xdb1(2:end) ... 
    - 20*Err_Ez_long_lg(2:end)  Xdb1(2:end) + 20*Err_Ez_long_lg(2:end) ]','color','r','LineWidth',2);


%% Screwed up JULIA. I was comapring the all data with CHAMP ap > 20
% now select JULIA only for ap > 20
load /nfs/satmag_work/mnair/projects/longp/aplist.mat;
load /nfs/satmag_work/mnair/projects/longp/alldays JULI_SEG ACE_SEG N_data TIME_SEG;
JULI_SEG = JULI_SEG.*24.366*1e-3; %mV/m
TIME_SEG = TIME_SEG + datenum(2000,1,1);
%% select only JULIA data for which ap >= 20 (This scripts also was used to make JULAI
% spectra for Ap < = 20

st = 1;
en = 72;
st_n = 1;
en_n = 72;
New_N_data = 0;
clear TIME_SEG_NEW JULI_SEG_NEW ACE_SEG_NEW
for i = 1:N_data,
    
    
 L = fday_ap >= TIME_SEG(st) & fday_ap <=  TIME_SEG(en);
                    mean_ap = mean(ap(L));
                    
                    if mean_ap >= 20,
                        fprintf('%10.5f %d %d %d %d\n',mean_ap, i, st, en, en-st);
                        TIME_SEG_NEW(st_n:en_n) = TIME_SEG(st:en);
                        JULI_SEG_NEW(st_n:en_n) = JULI_SEG(st:en);
                        ACE_SEG_NEW(st_n:en_n) = ACE_SEG(st:en);
                        st_n = st_n + 72;
                        en_n = en_n + 72;
                        New_N_data = New_N_data + 1;
                    end;
                    
                    st = st + 72;
                    en = en + 72;
                   
end;

%[Pxx_julia_apGT20,F_julia] = pwelch(JULI_SEG_NEW,hanning(72),0,72,1/(5*60));
[Pxx_julia_apLT20,F_julia] = pwelch(JULI_SEG_NEW,hanning(72),0,72,1/(5*60));
%% calculate the TF etc
L = ACE_SEG_NEW < 0;
temp = abs(ACE_SEG_NEW);
Eeff = 8 * temp ./ sqrt (8^2 + temp.^2);
ACE_SEG_NEW = Eeff;
ACE_SEG_NEW(L) = ACE_SEG_NEW(L) * -1;


[Txy_short,F_short] = tfestimate(ACE_SEG_NEW,JULI_SEG_NEW,hanning(72),0,72,1/(5*60)); %Txy is agains calculated just to get
[Cxy_short,F_short] = mscohere(ACE_SEG_NEW,JULI_SEG_NEW,hanning(72),0,72,1/(5*60));
New_N_data = length(ACE_SEG_NEW) / 72;
Err_short = sqrt( 1/(2*(New_N_data-1)) .* ( (1-Cxy_short)./Cxy_short ) ) .* abs(Txy_short);
Err_short_lg = 0.434 * Err_short ./ abs (Txy_short);

%plot the data

 errorbar(log10((1./F_short(2:end))./(3600)), log10(abs(Txy_short(2:end))),...
     (Err_short_lg(2:end)),(Err_short_lg(2:end)),'r', 'linewidth',2);
    
%% Re plotting ACE spectra
% The script in this cell was made to plot ACE spectra. First, the gaps in
% the data > 12 minutes are determined. 

load /nfs/satmag_work/mnair/projects/longp/aplist.mat;
load /nfs/satmag_work/mnair/projects/ace_tensor/acedata/ace_2000_2010.mat ace_all;

% finding gaps in ACE data with gaps > 12 minutes

nd = 1;
np = 0;

L = isnan(ace_all(:,2)) | isnan(ace_all(:,3));

for i = 1: length(ace_all) - 1,
    
    
    if  L(i) && L(i + 1)
        
        np = np + 1;
        
    else
        if ( ace_all(i,1) - ace_all(i - np ,1) > 12/1440 )
        % Mar 2012. One of the issues to be looked at is how these
        % gaps are affecting the spectra. A future improvement
        % should be to break the ACE data into pieces (just like
        % EEF data) and limit the Fourier analysis to within the
        % segements.
        data_index_ace(nd,2) = i;
        data_index_ace(nd,1) = i - np;
        nd = nd + 1;
        
    end;
    np = 0;
    
    end;
end;



L = isnan(ace_all(:,2));
y = interp1(ace_all(~L,1),ace_all(~L,2),ace_all(:,1),'spline');
ace_all(:,2) = y;


% filling long gaps with linear interpolation


for i = 1:length(data_index_ace),
    
%         ace_all(data_index_ace(i,1):data_index_ace(i,2),2) = 0;
%         ace_all(data_index_ace(i,1):data_index_ace(i,2),3) = 0;
%     
%     linear interpolation for the gaps > 12 minutes. Note that for gaps
%     less than 12 minutes, we use the defined interpolation method (such as
%     spline)
    ace_all(data_index_ace(i,1):data_index_ace(i,2),2) = ...
        interp1(ace_all([data_index_ace(i,1)-1 data_index_ace(i,2)+1],1), ...
        ace_all([data_index_ace(i,1)-1 data_index_ace(i,2)+1],2), ...
        ace_all(data_index_ace(i,1):data_index_ace(i,2),1));
    
    ace_all(data_index_ace(i,1):data_index_ace(i,2),3) = ...
        interp1(ace_all([data_index_ace(i,1)-1 data_index_ace(i,2)+1],1), ...
        ace_all([data_index_ace(i,1)-1 data_index_ace(i,2)+1],3), ...
        ace_all(data_index_ace(i,1):data_index_ace(i,2),1));
    
    
    
end;




len = 6144;

st = 1;
en = len;
st_n = 1;
en_n = len;
New_N_data = 0;
N_data = length(ace_all) / len;
clear TIME_SEG_NEW ACE_SEG_NEW
for i = 1:N_data*2-1
        
 L = fday_ap >= ace_all(st,1) & fday_ap <=  ace_all(en,1);
                    mean_ap = mean(ap(L));
                    
                    if mean_ap >= 20,
                        %fprintf('%10.5f %d %d %d %d\n',mean_ap, i, st, en, en-st);
                        TIME_SEG_NEW(st_n:en_n) = ace_all(st:en,1);
     
                        ACE_SEG_NEW(st_n:en_n) = ace_all(st:en,2);
                        st_n = st_n + len;
                        en_n = en_n + len;
                        New_N_data = New_N_data + 1;
                    end;
                    
                    st = st + len/2;
                    en = en + len/2;
                   
end;

[Pxx_ace,F_ace] = pwelch(ACE_SEG_NEW,hanning(len),0,len,1/(5*60));
loglog((1./(3600*F_ace(2:end))),abs(Pxx_ace(2:end)),'b-','LineWidth',2);

% Estimating spectrum using the new spectrum methods
% H = SPECTRUM.WELCH(WINNAME,SEGMENTLENGTH,OVERLAPPERCENT) 
% Hs=spectrum.welch('Hann',len,0);
% Hpsd = psd(Hs,ACE_SEG_NEW,'Fs',1/(5*60),'ConfLevel',0.95);

% The following script restrict the data to data with gaps < 12 minutes.
% However, the results are not very different from considering all data
% with gaps filled by linear interpolation.

% Find data index for good data 

% data_index_yes = [[0; data_index_ace(:,1)] [ data_index_ace(:,2) ;length(ace_all)]];
% 
% len = 1080;
% 
% st = 1;
% en = len;
% st_n = 1;
% en_n = len;
% New_N_data = 0;
% N_data = length(ace_all) / len;
% clear TIME_SEG_NEW ACE_SEG_NEW
% for i = 1:length(data_index_yes),
%         
%     st = data_index_yes(i,1);
%     en = data_index_yes(i,2);
%     
%     if en - st >= len,
%         en = st + len - 1; % consider only the first len data
%                     L = fday_ap >= ace_all(st,1) & fday_ap <=  ace_all(en,1);
%                     mean_ap = mean(ap(L));
%                     
%                     if mean_ap > 0,
%                         %fprintf('%10.5f %d %d %d %d\n',mean_ap, i, st, en, en-st);
%                         TIME_SEG_NEW(st_n:en_n) = ace_all(st:en,1);
%                         ACE_SEG_NEW(st_n:en_n) = ace_all(st:en,2);
%                         st_n = st_n + len;
%                         en_n = en_n + len;
%                         New_N_data = New_N_data + 1;
%                     end;
%                     
%     end;
%                    
% end;
% 
% [Pxx_ace,F_ace] = pwelch(ACE_SEG_NEW,hanning(len),0,len,1/(5*60));
% loglog((1./(3600*F_ace(2:end))),abs(Pxx_ace(2:end)),'b*-','LineWidth',2);
% 
% Estimating spectrum using the new spectrum methods
% H = SPECTRUM.WELCH(WINNAME,SEGMENTLENGTH,OVERLAPPERCENT) 
% Hs=spectrum.welch('Hann',len,0);
% Hpsd = psd(Hs,ACE_SEG_NEW,'Fs',1/(5*60),'ConfLevel',0.95);

%% plot TF magnitude in normal log scale (not DB)

load  /nfs/satmag_work/mnair/projects/longp/JULIA_GT20_FINAL
load /nfs/satmag_work/mnair/projects/longp/CHAMP_LONG_PERIOD_GT20_FINAL

%figure(2);
loglog((1./(3600*F_short)),abs(Txy_short),'bs','MarkerSize',5);
hold on;
loglog((1./(3600*F_long)),abs(Txy_long),'rs','MarkerSize',5);
set(gca,'FontSize',16);

for i = 1: length(Txy_short);
    plot([1./(3600*F_short(i)) 1./(3600*F_short(i)) ], [abs(Txy_short(i)) - ... 
        Err_short(i) abs(Txy_short(i)) + Err_short(i)],'b','LineWidth',1);
    hold on;
end;
%
    
for i = 1: length(Txy_long);
    plot([1./(3600*F_long(i)) 1./(3600*F_long(i)) ], [abs(Txy_long(i)) - ... 
        Err_long(i) abs(Txy_long(i)) + Err_long(i)],'r','LineWidth',1);
    hold on;
end;
    

set(gca,'XTick',[6/60,10/60,20/60,30/60,1,2,4,6,10, 30]);
set(gca,'XTickLabels',[' 6';'10';'20';'30';' 1';' 2';' 4';' 6';'10';'30']);
grid on;
axis([0.1 40 1e-3 1e-1 ]);
 aa=axis;
 text(aa(1),aa(3),'minutes','FontSize',16);
 text(aa(2),aa(3),'hours','FontSize',16);
 
 text(aa(1),aa(3),'JULIA','color','b','FontSize',16);
 text(aa(2),aa(3),'CHAMP','color','r','FontSize',16);
 
 % plot filter fit on this
 
 

[h_n] = freqz(b_n,a_n,F_all,1/(5*60));
h = loglog((1./(3600*F_all)),abs(h_n),'b','LineWidth',2);

%% compare the power spectra of ief ey computed by different approaches like
%welch average periodogram , pmtm and robust fit 

% Compute the spectra for ACE Ey using the Welch's average periodogram as
% described above

% compute pmtm and average
for i = 1:10800:length(ace_all),
[P1(icount,:),F1] = pmtm(ace_all(i:i+10800-1,2), 2, 10800, 1/(60*5));
icount = icount + 1;
end;
P1(1,1)
P1(100,100)
loglog((1./(3600*F1)),mean(P1,1),'b-','LineWidth',1);


% robust fitting

periods = 1./(F1(2:50:end));
spectra = [];
for i = 1:length(periods),

x = [cos(2*pi*ace_all(:,1)/( periods(i)/24)) sin(2*pi*ace_all(:,1)/( periods(i)/24))];	 
 
spectra_model_rem(i,:) = robustfit(x, ace_all(:,2));

end;

%% plot TF magnitude and phase of ALL Day JULIA 

load /nfs/satmag_work/mnair/projects/longp/alldays JULI_SEG ACE_SEG;
JULI_SEG = JULI_SEG.*24.366*1e-3; %mV/m
[Txy_short,F_short] = tfestimate(ACE_SEG,JULI_SEG,hanning(72),0,72,1/(5*60)); %Txy is agains calculated just to get
[Cxy_short,F_short] = mscohere(ACE_SEG,JULI_SEG,hanning(72),0,72,1/(5*60));
N_data = length(ACE_SEG) / 72;
Err_short = sqrt( 1/(2*(N_data-1)) .* ( (1-Cxy_short)./Cxy_short ) ) .* abs(Txy_short);
Err_phase_short = atan(sqrt( 1/(2*(N_data-1)) .* ( (1-Cxy_short)./Cxy_short ) ) ) * 180/pi;

%magnitude

%figure(2);
subplot(211);
loglog((1./(3600*F_short)),abs(Txy_short),'ks','MarkerSize',5, 'LineWidth',1);
hold on;
set(gca,'FontSize',16);

for i = 1: length(Txy_short);
    plot([1./(3600*F_short(i)) 1./(3600*F_short(i)) ], [abs(Txy_short(i)) - ... 
        Err_short(i) abs(Txy_short(i)) + Err_short(i)],'k','LineWidth',1);
    hold on;
end;

set(gca,'XTick',[6/60,10/60,20/60,30/60,1,2,4,6,10]);
set(gca,'XTickLabels',[' 6';'10';'20';'30';' 1';' 2';' 4';' 6';'10']);
grid on;
axis([0.1 10 1e-3 1e-1 ]);
 aa=axis;
% text(aa(1),aa(3),'minutes','FontSize',16);
% text(aa(2),aa(3),'hours','FontSize',16);
 
 %text(aa(1),aa(3),'JULIA','color','b','FontSize',16);

  ylabel('Magnitude ');
  
  % Phase
  subplot(212);
phase_data = angle(Txy_short) * 180/pi;

semilogx((1./(3600*F_short)),phase_data,'ks','MarkerSize',5,'LineWidth',1);
hold on;
set(gca,'FontSize',16);

for i = 1: length(phase_data);
    plot([1./(3600*F_short(i)) 1./(3600*F_short(i)) ], [phase_data(i) - ... 
        Err_phase_short(i) phase_data(i) + Err_phase_short(i)],'k','LineWidth',1);
    hold on;
end;


set(gca,'XTick',[6/60,10/60,20/60,30/60,1,2,4,6,10]);
set(gca,'XTickLabels',[' 6';'10';'20';'30';' 1';' 2';' 4';' 6';'10']);
grid on;
axis([0.1 10 -180 180 ]);
 aa=axis;
 text(aa(1),aa(3),'minutes','FontSize',16);
 text(aa(2),aa(3),'hours','FontSize',16);
 
 text(aa(1),aa(3),'JULIA','color','b','FontSize',16);

  ylabel('Phase ');
  
  
saveas(gcf, '/data/backup/mnair/sw_paper/Mag_Phase_transfer_function');
saveas(gcf, '/data/backup/mnair/sw_paper/Mag_Phase_transfer_function.eps','eps');

%% plot the model application for space weather paper (f)





load /data/backup/mnair/sw_paper/model_application_day_data
load /data/backup/mnair/sw_paper/kelley_india_drift_2002_apr_17 x y;


%Need to first process the fejer climate for Indian side

% fid = fopen('/nfs/satmag_work/mnair/projects/longp/jic_isr_fejer.txt','wt');
% 
% for i = 1:length(x),
%     
%     
%     f107d = 185;
%     doy = 107;
%     
%     in_lt = x(i) + 5.5; % 5.5 for India
%     
%     if in_lt > 24,
%         in_lt = in_lt - 24;
%     end;
%     fprintf(fid,'%3d %6.3f %6.3f \n', doy, in_lt, f107d);
% end;

subplot(311);

plot((ief_fday-floor(ief_fday))*24,ief_ey,(ief_fday-floor(ief_fday))*24,ief_ez, 'LineWidth', 2);

%plotting only IMF Ey
plot((ief_fday-floor(ief_fday))*24,ief_ey, 'LineWidth', 2);

axis([0 24 -inf inf]);

xticks = 0:2:24 ;

set(gca,'XTick', xticks);


%     xticks_lt = xticks - 5;
%
%     xticks_lt (xticks_lt <= 0) = xticks_lt (xticks_lt <= 0) + 24;
%
%     set(gca,'XTickLabel',  reshape(sprintf('%2.0f',xticks_lt),[2,length(xticks_lt)])' );

set(gca,'FontSize',16);
xlabel('UT (Hours)');
ylabel('IEF mV/m');
legend('IEF Ey');
text(0,0,'ACE IEF Ey', 'FontSize',20)

title(['Date = ' datestr(A(i)) ', ' sprintf('Mean Ap = %7.0f', mean_ap)]);
axis([12,24,-inf,inf]);

subplot(312);

if jic_flag == 1 && length(jicamarca_eef_fday) > 10,

    % runmean on window 2*M + 1
    jic_data = runmean(jicamarca_eef_piece_old - jicamarca_quiet_day, 2) ;
    
    plot((jicamarca_eef_fday_old - floor(jicamarca_eef_fday_old ) ) * 24 + 0.16, ...
        jic_data , 'b','LineWidth', 2);
    hold on;
    
end;




plot((predicted_eef_ut_fday - floor(predicted_eef_ut_fday))*24, ...
    jicamarca_eef_lt .* predicted_eef_ut,'r', 'LineWidth', 2);


axis([0 24 -inf inf]);

xticks = 0:2:24 ;

set(gca,'XTick', xticks);


xticks_lt = xticks - 5;

xticks_lt (xticks_lt <= 0) = xticks_lt (xticks_lt <= 0) + 24;

set(gca,'XTickLabel',  reshape(sprintf('%2.0f',xticks_lt),[2,length(xticks_lt)])' );

set(gca,'FontSize',16);
xlabel('Local Time (Hours)');
ylabel('EEF mV/m');
legend('Observed', 'Predicted');

text(13,0,'South America', 'FontSize',20)
hold off;
 axis([12,24,-inf,inf]);
subplot(313);

%Calculate the EEF in Indian sector following the eqution
% Drift  = 5.2889 + 0.1947DH + 0.0001DH^2- 0.0000021DH^3 by
%Anderson et al 2004 Space Weather and using B = 37000 nT (2004,
%trivandrum, IGRF at 150 km)

eef_india = (5.2889 + 0.1947*india_delh + 0.0001*(india_delh).^2 -...
    0.0000021*(india_delh).^3) * 37000 / 1e6;

%         plot((ind_fday - floor(ind_fday))*24 + phase_delay_correction(n_phase_delay_corr)/1440 ...
%         , eef_india - mean(eef_india), 'LineWidth', 2);
%plot((ind_fday - floor(ind_fday))*24, eef_india - mean(eef_india), 'LineWidth', 2);
%hold on;
ind_drift = load('/data/backup/mnair/sw_paper/ind_drift_fejer.out');
%plot(x,(y-58)*30*1e-3 + ind_drift(1:end-1,2)*30*1e-3 ,'k');
plot(x +0.16,(y-58)*30*1e-3 - mean((y-58)*30*1e-3),'b','LineWidth', 2);
hold on;
plot((predicted_eef_ut_fday - floor(predicted_eef_ut_fday))*24, india_eef_lt.*predicted_eef_ut,'r', 'LineWidth', 2);
%plot((predicted_eef_ut_fday - floor(predicted_eef_ut_fday))*24, -1*predicted_eef_ut,'r', 'LineWidth', 2);
axis([0 24 -inf inf]);

xticks = 0:2:24 ;

set(gca,'XTick', xticks);
% This is to increment the phase delay array index. This should be
% used for plotting a few data sets. For large data scrolling,
% remove this and its other references.
n_phase_delay_corr = n_phase_delay_corr + 1;

xticks_lt = xticks + 5.5;

xticks_lt (xticks_lt >= 24) = xticks_lt (xticks_lt >= 24) - 24;

set(gca,'XTickLabel',  reshape(sprintf('%4.1f',xticks_lt),[4,length(xticks_lt)])' );

set(gca,'FontSize',16);
xlabel('Local Time (Hours)');
ylabel('EEF mV/m');
legend('Observed','Predicted');
text(13,0,'India', 'FontSize',20)
 axis([12,24,-inf,inf]);

set(gcf,'Position', [670   389   709   701]);
saveas(gcf, '/data/backup/mnair/sw_paper/20020417_obs_pre.fig','fig');
saveas(gcf, '/data/backup/mnair/sw_paper/20020417_obs_pre.eps','eps');


%% Plot the ACE, JULIA and CHAMP EEF Spectra. T
% The idea is that this section should produce the figures 2 (or 2 and 3)
% of the CHAMP JGR paper 

% Load the ACE spectra when Ap is between 0 and 20
load /nfs/satmag_work/mnair/projects/longp/ACE_AP_0_20 Pxx_ace F_ace;
P_ace_ap_0_20 = Pxx_ace;

% Load the ACE spectra when Ap is between 20 and inf
load /nfs/satmag_work/mnair/projects/longp/ACE_AP_GT_20 Pxx_ace F_ace;
P_ace_ap_20_inf = Pxx_ace;

% load the JULIA data for AP 0-20 and 20-Inf
load /nfs/satmag_work/mnair/projects/longp/ACE_JULIA_SEG_AP_0_20_AND_GT_20 Pxx_julia_apLT20 Pxx_julia_apGT20 F_julia

% Load the CHAMP EEF data 0-20
%load /nfs/satmag_work/mnair/projects/longp/ACE_CHAMP_SEG_AP_0_20 Pxx F_long
load /nfs/satmag_work/mnair/projects/longp/ACE_CHAMP_SEG_AP_0_20_512 Pxx F_long
P_champ_ap_0_20 = Pxx;
% Load the CHAMP EEF data 20-inf
load /nfs/satmag_work/mnair/projects/longp/ACE_CHAMP_SEG_AP_20_1000_512 Pxx F_long
P_champ_ap_20_inf = Pxx;

% Now plot the stuff

% loglog((1./(3600*F_ace(2:end))),abs(P_ace_ap_0_20(2:end)),'k--','LineWidth',2);
% hold on;
% loglog((1./(3600*F_ace(2:end))),abs(P_ace_ap_20_inf(2:end)),'k-','LineWidth',2);
set(gca,'FontSize',16);
set(gcf, 'PaperPositionMode', 'auto');
grid on;
 set(gca,'YMinorGrid','off');
 set(gca,'XMinorGrid','off');
ylabel('Power Spectral Density $[\frac{mV^2}{m^2}] / Hz$', 'Interpreter', 'Latex');
xlabel('Period (Hours)', 'Interpreter', 'Latex');

legend('Ap < 20','Ap > 20','Interpreter', 'Latex');
title('IEF E_y');
saveas(gcf, '/nfs/satmag_work/mnair/projects/longp/Figures/ACE_Spectra','fig');
saveas(gcf, '/nfs/satmag_work/mnair/projects/longp/Figures/ACE_Spectra','png');
saveas(gcf, '/nfs/satmag_work/mnair/projects/longp/Figures/ACE_Spectra.eps','psc2');

% plot the EEF from CHAMP and JULIA

loglog((1./(3600*F_julia(2:end))),abs(Pxx_julia_apLT20(2:end)),'r--','LineWidth',2);
hold on;
loglog((1./(3600*F_julia(2:end))),abs(Pxx_julia_apGT20(2:end)),'b--','LineWidth',2);

%only plot data for periods > 6 hours (to remove the aliasing effects)
%and 100 hours
L = (1./(3600*F_long(1:end))) > 6.3 & (1./(3600*F_long(1:end))) < 700;

loglog((1./(3600*F_long(L))),abs(P_champ_ap_0_20(L)),'r','LineWidth',2);

loglog((1./(3600*F_long(L))),abs(P_champ_ap_20_inf(L)),'b','LineWidth',2);

set(gca,'FontSize',16);
set(gcf, 'PaperPositionMode', 'auto');
grid;
 set(gca,'YMinorGrid','off');
 set(gca,'XMinorGrid','off');
ylabel('Power Spectral Density $[\frac{mV^2}{m^2}] / Hz$', 'Interpreter', 'Latex');
xlabel('Period (Hours)', 'Interpreter', 'Latex');
legend('JULIA Ap < 20','JULIA Ap > 20','CHAMP Ap < 20','CHAMP Ap > 20');
title('Ionospheric Zonal Electric Field');
saveas(gcf, '/nfs/satmag_work/mnair/projects/longp/Figures/EEF_Spectra','fig');
saveas(gcf, '/nfs/satmag_work/mnair/projects/longp/Figures/EEF_Spectra','png');
saveas(gcf, '/nfs/satmag_work/mnair/projects/longp/Figures/EEF_Spectra.eps','psc2');
%print -deps2 -loose '/nfs/satmag_work/mnair/projects/longp/Figures/EEF_Spectra.ps'
